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INTRODUCTION 


The  idea  of  estimating  the  state  of  a  system  based  on  observed  data  has 

been  the  focal  point  of  much  research.  The  concept  of  least  squares  and  curve 

fitting  was  introduced  by  Gauss1  in  the  early  1800s.  In  more  recent  times, 

2 

Norbert  Wiener  was  asked  to  solve  an  important  specific  World  War  II  problem: 
"How  can  optimum  properties  for  servomechanisms  be  formulated  and  specified?" 
Wiener  formulated  a  general  solution  to  this  problem  based  on  rigorous  proba¬ 
bilistic  approaches.  Among  his  works,  "The  Wiener  Filter"  is  widely  acclaimed 
today  as  a  cornerstone  of  modern  estimation  theory.  Over  the  past  decade  and 

a  half,  modern  sequential  estimation  theory  has  been  coming  of  age.  Kalman3 

4 

and  Kalman  and  Bucy  brought  modern  estimation  theory  into  the  limelight  with 
their  historical  publications,  "A  New  Approach  to  Linear  Filtering  and  Pre¬ 
diction  Theory"  and  "New  Results  in  Linear  Filtering  and  Prediction  Theory." 
Scientists,  engineers,  statisticians,  mathematicians,  and  technical  personnel 
faced  with  the  problem  of  estimating  the  behavior  or  monitoring  the  state 
variables  of  a  physical  system  find  themselves  relying  heavily  on  modern 
estimation  theory.  The  problem  is  not  an  may  one  in  that  the  choice  of  the 
estimation  technique,  in  many  cases,  is  problem-dependent.  There  are  many 
criteria  available  to  specify  the  methodology  for  estimating  a  parameter  or 

variable.  Ho  and  Lee  discuss  three  choices:  a)  most  probable  estimate;  b) 

£ 

conditional  mean  estimate;  and  c)  minimax  estimate.  Jacquot  presents  a 

1.  Gauss,  K.  F.,  Theory  of  Motion  of  the  Heavenly  Bodies,  Dover,  1963. 

2.  Wiener,  N.,  The  Extrapolation,  Interpolation,  and  Smoothing  of  Stationary 
Time  Series  with  Engineering  Applications,  MIT  Press,  1949. 

3.  Kalman,  R.  E.,  A  New  Approach  to  Linear  Filtering  and  Prediction  Theory, 
Trans  ASME,  (Journal  of  Basic  Engineering)  ser  D,  vol  82,  pp  35-45,  March 
1960. 

4.  Kalman,  R.  B. ,  and  R.  S.  Bucy,  New  Results  in  Linear  Filtering  and  Pre¬ 
diction  Theory,  Trans  ASME,  (Journal  of  Basic  Engineering) ,  ser  D,  vol 
83,  pp  95-108,  March  1961. 

5.  Ho,  V.  C.,  and  R.  C.  K.  Lee,  A  Bayesian  Approach  to  Problems  in  Sto¬ 
chastic  Estimation  and  Control,  IEEE  Trans  Automatic  Control,  vol  AC-9, 
pp  333-339,  October  1969. 

6.  Jacquot,  R.  G.,  Modern  Digital  Control,  Marcel  Dekker,  Inc,  1981. 
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detailed  derivation  of  the  classical  discrete-time  vector  Kalman  filter  where 
the  system  of  interest  is  governed  by  the  stochastic  matrix-vector  difference 
equation 


x(k  +  1)  *  Ax( k)  +  w(k)  +  Bu(k)  (1) 

with  the  measurement  process  defined  by 

z(k)  =  Hx(k)  +  v(k).  (2) 

In  equation  (1),  x(k)  is  an  n  X  1  vector;  A  is  an  n  X  n  state  transition 
matrix;  B  is  an  n  X  j  distribution  matrix;  and  w(k)  and  u(k)  are  n  X  1  and  j  X 
1  vectors,  respectively.  The  stochastic  sequence  w(k)  is  defined  as  the  state 
or  process  noise.  In  equation  (2),  the  measurement  equation,  z(k)  is  an  r  X  1 
vector  of  observations,  H  is  an  r  X  n  measurement  matrix,  and  v(k)  is  the  r  X 
1  vector  of  measurement  noise. 

The  assumptions  will  be  made  that  the  noises  w(k)  and  v(k)  are  sta¬ 
tionary,  Gaussian  random  sequences.  The  problem  is  to  obtain  the  best  esti¬ 
mate  of  x(k)  such  that  the  conditional  mean  square  error 

J  =  1/2  E[xT(k)  x(k)|Zk]  (3) 


is  minimized.  The  estimation  error  x(k)  is  given  as 

x(k)  »  x(k)  -  x(k).  (4) 

The  notation  Z*  =  {z(1),  z(2),  ...,  z(k)}  is  used  to  represent  the  past 
observation  or  measurement  set.  The  vector  x(k)  is  the  estimate  of  the  true 
vector  x(k).  The  conditional  mean  square  error  of  (3)  is  conditional  on  the 
past  observations  or  measurements.  The  estimation  process  is  formalized  in 
the  following  algorithm: 

x(k)  *  Ax( k  -  1)  +  q  +  Bu(k  -  1)  (5) 
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P(k)  -  A$(k  -  1)AT  +  Q 


(6 


Observation  Residual: 

y(k)  -  *(k)  -  Hx(k) .  (7) 

Kalman  Gain: 

K(k)  -  P(k)  HT[HP(k)HT  +  R]"1.  (8) 

Estimation  Equations: 

Measurement  update- 

x(k)  -  x(k)  +  K(k)[y(k)  -  r  -  HBu( k  -  1)  -  Hq]  (9) 

Covariance  update— 

lP(k)  -  P(k)  -  K( k)HP(k)  •  (10) 

A  sequence  which  is  of  extreme  importance  in  the  field  of  estimation  is 
termed  the  innovations  sequence  and  is  defined  as  the  observation  residual 
given  by  equation  ( 7 ) . 

This  development  of  the  estimation  problem  is  the  classical  problem  for 
known  statistics:  ie ,  R  and  Q  are  known  quantities:  also,  the  noise  sequences 
have  the  properties 

tlv(k)]  -  r  (11) 


and 


E(w(k)l  -  q 


( 12 


The  classical  problem  also  assumes  that  if  the  system  is  driven  by  a  deter-* 
ministic  forcing  function,  it  is  completely  known  for  all  time. 


The  noise  covariance  matrix  Q  is  defined  by 

Q  ■  E  |[w(k)  -  q]  [w(k)  -  q]T|  M3) 

where  w(k)  is  the  random  forcing  sequence  of  expression  (1)  and  q  as  defined 
in  relation  (12),  wnere  the  mean  has  been  assumed  not  to  be  a  function  of  time 
or  k.  The  observation  noise  covariance  matrix  R  is  given  as 

R  -  E  -j[v(k)  -  r]  (v(k)  -  r]T  }  (14) 

where  v(k)  is  the  measurement  noise  sequence  and  r  is  the  expected  value  of 
v(k). 


Figure  1  presents  a  block  diagram  of  a  system  model,  measurement  system, 
and  discrete  Kalman  filter.  The  system  model  is  a  discrete  representation  of 
a  continuous  system  being  observed  at  discrete  times  by  a  measurement  system. 
The  assumptions  are  that  u(k)  and  w(k)  are  constant  over  the  sampling  inter¬ 
val  ;  ie , 

w(t)  *  w(k);  kT<_  t  <.(k  +  1)T  (15) 


and 


u(t)  -  u(k) ;  kT<_  t  <(k  +  1)T 
where  T  is  the  sampling  interval. 
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Figure  1 .  System  model  and  conventional  discrete  Kalman  filter. 


A  further  assumption  is  that  knowledge  of  the  forcing  function  exists;  thus 
the  Kalman  filter,  as  illustrated  in  figure  1,  reflects  this  knowledge.  Con¬ 
sider  now  the  case  where  a  system  is  being  forced  by  forcing  functions  u(k) 
and  w(k)  and  the  noise  corrupted  output  is  observed  by  a  sensor.  The  deter¬ 
ministic  forcing  function,  u(k),  is  fast  varying  with  respect  to  time.  It  is 
desired  to  formulate  estimates  of  the  value  of  the  state  variables  in  a  timely 
manner  by  using  a  Kalman  filter;  however,  for  this  case,  knowledge  of  the 
system  forcing  function,  u(k),  is  unknown  to  the  filter.  A  problem  of  equal 
importance  is  the  case  where  the  statistics  of  the  process  noise,  w(k),  are 
unknown.  An  excellent  treatment  of  estimation  in  the  presence  of  unknown 
noise  statistics  is  presented  in  Myers  and  Tapley  .  Empirical  estimators  are 
developed  in  reference  [7]  which  estimate  the  noise  statistics. 


For  the  case  at  hand,  state  estimation  without  knowledge  of  the  deter¬ 
ministic  forcing  functions,  several  modifications  are  incorporated  into  the 
estimation  process.  These  include  adaptive  weighting  of  the  elements  of  the 


7.  Myers,  K.  A.  and  B.  D.  Tapley,  Adaptive  Sequential  Estimation  with 
Unknown  Noise  Statistics,  IEEE  Trans  Automatic  Control,  vol  AC-21, 
pp  520-523,  August  1976. 


conventional  Kalman  gain  and  covariance  matrices,  as  well  as  robust  statis¬ 
tical  smoothing  of  the  estimates  made  by  the  adaptive  Kalman  filter  using  the 
modified  gain  and  covariance  matrices.  Figure  2  illustrates  the  estimation 
process  with  the  modifications  incorporated.  The  assumptions  are  that  the 
noises  w(k)  and  v(k)  are  stationary,  Gaussian,  random  sequences.  The  modified 
gain  matrix  A[K(k)]  is  defined  in  the  next  section  as  well  as  the  robust 
statistical  smoothing  procedure. 


Figure  2.  Adaptive  robust  estimation. 


Two  concepts  are  investigated  in  this  paper;  1)  adaptive  weight  functions 
for  the  Kalman  filter  gain  and  error  covariance  matrices,  and  2)  robust 
smoothing  of  the  estimated  state  variables. 

Robust  smoothing  can  be  applied  to  both  the  observed  and  nonobserved 
state  variables.  The  concept  of  robust  statistical  smoothing  of  nonobserved 
variables  is  addressed  in  Groutage10,  where  it  is  applied  to  the  maneuvering 
target  tracking  problem. 

Basically,  two  robust  statistical  procedures  are  used  for  the  robust 
smoothing  of  the  estimated  state  variables: 

1)  robust  measure  of  spread,  and 
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^  •- 


i. 


2)  robust  measure  of  location* 

8  9 

Details  on  these  robust  statistics  are  found  in  Bickel  and  Huber  .  The 
concepts  of  using  adaptive  weights  for  Kalman  filter  gain  and  covariance 
matrices  are  presented  in  [10]. 

ROBUST  ESTIMATION  OF  OBSERVED  STATE  VARIABLES  BY  USING  ADAPTIVE  WEIGHTS 
FOR  GAIN  AND  ERROR  COVARIANCE  MATRICES 


The  concepts  presented  here  are  empirical  in  nature;  ie,  they  are  based 
on  observations  and  experimental  data.  An  intuitive  perception  of  these  con¬ 
cepts  can  be  obtained  by  examining  the  Kalman  filtering  algorithm  for  a  scalar 
parameter .  Consider  the  scalar  system  with  process  noise  w(k) 

x(k  +  1)  -  a  x(k)  +  w(k)  ( 17) 

and  the  associated  measurements 

z(k)  *  x(k)  +  v(k) .  (18) 

The  covariance  associated  with  the  process  noise  is 

E[w2(k)]  -  Q  (19) 

and  the  measurement  noise  covariance  is 

E(v2(k)]  «  R.  (20) 


8.  Bickel,  P.  J.,  One-Step  Huber  Estimates  in  the  Linear  Model,  Journal  of 
the  American  Statistical  Association,  vol  70,  no  350,  pp  428-433,.  June 
1975. 

9.  Huber,  P.  J.,  Robust  Estimation  of  a  Location  Parameter,  The  Annals  of 
MathesMtical  Statistics,  vol  35,  no  1,-  pp  73-101,  March  1964. 

10.  Groutage,  F.  D.,  A»,  ve  oust  Sequential  Estimation  with  Application 
to  Tracking  a  Maneuv-.ng  Target,  PhD  Dissertation,  Department  of  Elec¬ 
trical  Engineering,  university  of  Wyoming,  May  1982. 
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The  noises  w(k)  and  v(k)  are  zero  mean,  Gaussian  random  sequences.  The  Kalman 
filter  equations  for  this  scalar  case,  which  formulate  estimates  of  the  single 
state  variable  x(k),  are: 


Propagation  Equations 
x(k)  =  a  x(k  -  1) 

P(k)  =  a2$(k  -  1)  +  Q 
Gain  Relationship 


c(k) 


P(k) 

P(k)  +  R 


(21) 

(22) 


(23) 


Estimation  Equations 

x(k)  =  x(k)  +  c(k)[z(k)  -  x(k)]  (24) 

£(k)  -  [1  -  c(k)]P(k).  (25) 

Equation  (26)  for  estimating  the  state  variable  x(k)  contains  interesting 
information  concerning  the  estimation  process.  Note  that  the  scalar  gain, 
c(k),  is  bounded  from  above  and  below  as 

0  <  c(k)  <  1.  (26) 

For  the  case  when  c(k)  =  0,  equation  (24)  indicates  that  total  faith  is  placed 
in  the  estimation  process.  In  fact,  the  measurements  are  ignored  and  the  pre¬ 
vious  estimate  is  the  updated  estimate.  Now  consider  the  case  where  c(k)  =  1. 
For  this  limit,  equation  (24)  indicates  that  there  will  be  no  faith  in  the 
estimation  process.  In  fact,  the  current  measurement  for  the  updated  estimate 
of  the  state  variable  is  used,  with  these  concepts  in  mind,  the  idea  of  a 
pseudogain,  a(k),  is  investigated.  Let 


a(  k) 


1  + 


-8(k) 

e  c(k) 


-v8(k) 

e 


(27) 


where  v  and  0(k)  are  parameters  to  be  established  and  v  can  be  a  constant. 
The  observation  residual,  y(k),  is  defined  as  the  difference  between  the 
measurement  and  the  propagated  state  estimate 


y(k)  -  z(k)  -  x(k) .  ( 2£ 

Now  consider  the  recursive  sample  mean  of  the  observation  residual  sequence 
and  the  recursive  sample  variance  of  the  residual  sequence.  Let  the  sample 
space  be  N  and  let  y(k)  designate  the  recursive  sample  mean  of  the  observa- 
tion  residual  sequence;  ie. 


.  * 

yoo  ■  jj  Y 


i  j=k-Nt+  1 


Let  0^  (k)  designate  the  recursive  sample  variance  (see  appendix  A)  of  the 
observation  residual  sequence;  ie. 


o  2(k)  "  3  2<k  -  D  +  „  1  ,  Uy(k)  ~  y(k)]2  -  [y(k  -  N)  -  y(k) ] 2 

y  y  ■  I  * 

+  [y(k)  -  y(k  -  N)]2}. 


If  the  parameter  B(k)  of  equation  (27)  is  chosen  in  the  following  manner 

B(k)  -  y®y2(k)  (31) 

then  a(k)  is  a  function  of  the  dispersion  of  the  residuals  y(k).  The  value  of 
the  constant,  y  (a  weighting  factor),  must  be  determined.  Equation  (27)  is 
now  interpreted.  Note  that  in  the  limit  as  0(k)  •*  0,  the  pseudogain  ap¬ 
proaches  the  Kalman  filter  gain  or 


lim 

8(k)*0 


a(k)  -  c(k) 
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and  for  the  limit  as  8(k)  ♦  • 


lim  a ( k)  =  1. 
8(k)— 


Thus  for  a  small  dispersion  of  the  residuals,  the  gain  is  the  optimum  Kalman 
gain  c(k).  For  a  large  dispersion  of  the  residuals,  no  faith  is  given  to  the 
estimation  process. 

From  practical  experiments ,  it  was  found  that  a  robust  weighting  function 
for  the  propagated  error  covariance  P(k)  was  also  required.  The  modified 
error  covariance  was  defined  as 

0 • < k)  -  P(k)  +  f • [ 1  -  a { k) ]  (34) 


where 


P'(k)  *  o(k)P(k) . 


For  the  scalar  case,  where  o(k)  is  defined  as  e 


A  2 

-vl0y  (k) 


-vo  2(k)  o1o  2 ( k )  -vo  2(k) 

0 ’ ( k)  -  e  y  P(k)  +  e  y  [1  -  e  y  J. 


^  2..  . 

a,o  (k) 

The  quantity,  e  1  is  limited  to  some  a  priori  upper  bound,  T.  This  is 

illustrated  in  figure  3.  Note  that  for  a  small  dispersion  of  the  residuals 
2 

(3  (k)  approaches  zero],  the  modified  error  covariance  0'(k)  of  expression 

(38)  approaches  the  Kalman  propagated  error  covariance  P(k)  or 


0'(k)  -  P(k) 


3  2(k)  *  0 

y 
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Figure  3.  Influence  curve  of  sample  variance  for  observation  residual. 


For  the  other  limit,  that  is,  when  the  dispersion  of  the  residuals  is  large 
[3^200  becomes  large],  the  term  8 '  (k)  approaches  the  a  priori  upper  bound,  or 

lim  e'(k)  -  T  (38) 

»  2 

o  (k)  ♦  • 

y 

Note  that  the  gain,  c(k),  is  also  bounded  when  the  dispersion  of  the  residuals 
is  large;  ie. 


u"  c(l°  ■  tTT 

3  Z(k)  ♦  - 

y 


(39) 


These  concepts  for  the  scalar  case,  robust  weighting  of  the  Kalman  gain  and 
error  covariance  matrices,  are  expanded  to  the  vector  case.  Note  that  for  the 
scalar  case,  the  gain  and  covariance  matrices  are  also  scalars. 


ADAPTIVE  GAIN  MATRIX  WEIGHTING 

To  be  useful  in  an  adaptive  filter ,  expression  (27)  must  be  expanded  to 
the  vector  dynamic  case*  This  is  accomplished  with  the  A[K(k)]  matrix  (for 
linear  measurements  or  pseudolinear  measurements)  similar  to  the  a(k)  function 
for  the  scalar  case  of  equation  (27)*  When  the  measurements  are  linear, 

A [K(k) ]  will  replace  the  Kalman  gain  matrix  K(k)  in  the  estimation  algorithm. 

Consider  a  matrix  A[K(k)]  defined  as 

A [K(k) ]  -  [I*  +  EK(k)  -  F]  (40) 

where  I',  E  and  F  are  n  X  r,  n  X  n,  and  n  X  r  matrices,  respectively.  Recall 
that  x(k)  and  z(k)  are  n  X  1  and  r  X  1  vectors,  respectively.  The  individual 
matrices  of  equation  (42)  can  best  be  explained  through  example.  Consider  the 
linear  system  of  figure  4.  This  system  is  driven  by  a  deterministic  forcing 
function,  u(t),  and  process  noise,  w(t).  The  output  is  observed  discretely 
with  a  sensor  system.  The  measurements,  z(k),  made  by  the  sensor  are  cor¬ 
rupted  by  a  noise  sequence  v(k).  The  measurement  equation  is  thus 

z(k)  -  x^k)  +  v(k).  (41) 


f SYSTEM  1  p 

|  FORCING  wit  I  |  . 


Figure  4.  Linear  system  example. 
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In  matrix  fora. 


s(k) 


For  this  axaaple,  the  I1,  E,  and  F  matrices  are  defined  as  follow: 


(42) 


(43) 


(44) 


(45) 


where  0(k)  is  as  defined  in  relation  (31).  It  should  also  be  recalled  that 
for  this  case  the  K(k)  matrix  is  a  2  X  1  vector. 


ADAPTIVE  ERROR  COVARIANCE  MATRIX  WEIGHTING 


The  equivalent  matrix  formulation  of  expression  (36)  is 

0’(k)  -  P'(k)  +  F'(I  -  I).  (46) 

Again,  the  individual  matrices  of  (46)  can  best  be  explained  through  applica¬ 
tion  to  the  linear  problem  of  figure  4.  The  propagated  error  covariance 
aiatrix,  P(k),  is  partitioned  into  three  separate  matrices.  The  following 
notation  is  introduced 

P(k)  -  L(k)  ♦  D(k)  ♦  V(k)  (47) 

where  L(k)  is  a  lower  triangular  matrix  with  zeros  on  the  diagonal,  D(k)  is  a 
diagonal  matrix,  and  V(k)  is  an  upper  triangular  matrix  with  zeros  on  the 
diagonal.  The  P'(k)  matrix  is  defined  as 
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& 

I 


p'(lc)  -  V(k)  +  I  D(k)  +  L(k) 


(48) 


where  the  £  matrix  is  a  diagonal  matrix.  The  individual  elements  of  the  £ 

a  2 

"Vi°y  <*>  O 

matrix  are  of  the  form  e  Ji  where  o  is  a  sample  variance  of  the  ith 

yi 

observation  residual  sequence.  For  the  case  when  all  state  variables  are 
measured,  then  all  of  the  diagonal  elements  of  £  would  be  of  the  form 

2, 

K) 

For  the 

problem  at  hand,  where  only  discrete  measurements  of  the  output  variable, 

Xj (k) ,  are  available 


’Vy  2(k)  “Vv  "(k) 

e  i  .  The  ith  diagonal  element,  £^,  is  e  i 


£  - 


a  2 ... 


0 

1 


(49) 


and 


(50) 


For  a  case  where  measurements  of  both  state  variables,  x^k)  and  x2(k), 
were  available  the  £  and  F'  matrices  would  be  as  follows: 


o 

e  2, 


-vo  (k) 
e  2  y2 


(51) 


F*  - 


0 


(52) 


A  2 , 

a  o  (k) 
0  •  y2 


Vy  2(k)  «23  2(k) 

The  quantities  e  1  and  e  2  would  be  limited  similarly  as  shown 
in  figure  3  to  upper  bounds  of  T^  and  T2,  respectively. 

Both  the  adaptive  Kalman  gain  and  error  covariance  weighting  procedures 
were  incorporated  into  the  algorithm  of  equations  (5)  through  (10).  It  was 
determined  through  experimentation  that  these  robust  adaptive  procedures  were 
most  effective  if  they  were  activated  only  after  the  sample  variance  of  the 
innovations  sequence  reached  a  predetermined  threshold  level.  As  a  rule  of 
thumb,  the  threshold  was  taken  to  be  one  and  one- half  the  anticipated  value  of 
the  standard  deviation  of  the  observation  noise. 


ROBUST  SMOOTHING 


The  estimates  of  the  state  variables  made  by  the  Kalman  filter  with  the 
modified  gain  and  covariance  matrices  contain  periodic  outliers.  This  is  a 
result  of  the  sampling  procedure  and  the  way  in  which  the  sample  statistics  of 
the  innovations  sequence  are  utilized  to  formulate  weights  for  the  elements  of 
the  gain  and  covariance  matrices.  To  alleviate  the  outliers  in  the  state 
estimates,  a  robust  statistical  smoothing  procedure  was  incorporated  into  the 
estimation  procedure.  The  robust  smoother  uses  a  regression  procedure  in  the 
following  manner.  Consider  n  samples  of  estimates  of  the  ith  state  variable 
x^(k)  where  the  samples  are  defined  as  the  set 

X  -  {^(k  -  n-1),  ^(k  -  n),  ^(k  -  n+1),  ...x^k)).  (53) 

It  is  desired  to  find  a  weighted- least-squares  solution  for  the  straight-line 
regression  fit  through  the  n  samples  of  the  estimates  of  the  ith  state  vari¬ 
able,  x^,  over  the  discrete  interval  spanned  from  discrete  time  k-n-1  to 
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time  k.  The  specific  weighted-least-squares  solution  for  the  straight-line 
regression  case  (ie,  Y  *  0Q  +  61  K  +  E)  is  given  by  the  formulas 


16 


where 


end  is  the  ith  observation  with  %  the  estimate  of  location  based  on  n 
observations , 


£  [V(e  )t  ) 
i-1 
M 

2><o 

1  * 


A  robust  measure  of  scale  is  dsfined  in  [12]  as 

(Interquartile  Distance) 

*  2(0.6745)  11 

where  the  interquartile  range  is  dsfined  as  the  third  quartile  minus  the  first 
quartile  and  thus  gives  the  length  of  the  interval  in  which  the  middle  50 
percent  of  the  data  fall.  For  samples  thet  arise  from  Gaussian  distributions. 


s  is  an  estimate  of  o,  the  standard  deviation.  The  value  of  the  constant  c  is 
arbitrary.  To  have  a  feel  for  the  range  of  the  value  of  c,  note  that  with 

c  -  4 

cs  -  4o.  (62) 

Discrete  valuee  of  the  weighted- least-squares  solution  along  the  regres¬ 
sion  line  are  obtained  from  the  relationship 


t(R)  -  |Q  ♦  e,  R 


for  Jc  ■  1,  2,  . . . ,  n 


12.  Launer ,  R.  L. ,  and  G.  N.  Wilkinson,  Robustness  in  Statistics,  Academic 
Frees,  New  York,  1979. 


where  8^  and  81  are  defined  by  equations  (55)  and  (54).  The  vertical  dis¬ 
tances  from  regression  line  of  (63)  to  individual  data  points  at  the  n  dis¬ 
crete  tines  are  called  residuals  and  defined  by  the  relationship 


r(j)  -  (k  -  j  -  1)  -  ?(j) 


for  j  ■  1,  2,  . , 


The  above  expressions  of  (63)  and  (64)  are  used  in  formulating  the  robust 
smoothing  procedure  illustrated  in  figure  2.  This  robust  procedure  is  imple¬ 
mented  by  using  the  weighted- least-squares  solution  of  (63)  to  project  n  -  1 
past  values  of  the  estimates  (as  formulated  by  the  adaptive  filter)  of  the  ith 
state  variable  up  to  the  present  discrete  time,  t  -  k.  The  n  -  1  past  values 


of  the  estimates  of  the  ith  state  variable;  ie. 


{xt(k  -  n  -  1),  x1  (k  -  n),  x^k  -  1)>  (65) 

are  projected  to  discrete  time  t  ■  k  and  define  n  values  of  the  random  vari¬ 
able,  x j ( k) j  ie, 

x^ ( k)  -  ?(k)  +  r( j)  for  j  -  t,  2,  ...,  n.  (66) 

The  newly  formed  random  variable,  x^(k),  is  smoothed  by  using  the  relationship 


W.x. (k) 


*<k) 


izi _ 

n  u 

£  "j 

j-i 


where  ft^'fk)  is  the  smoothed  value  of  the  estimated  value  of  the  ith  state 
variable  as  generated  by  the  modified  gain  and  covariance  Kalman  filter.  A 
new  estimate,  at  discrete  time  k+1,  of  the  ith  state  variable  is  generated, 
x^'(k+1),  which  is  subsequently  ssxx>thed  by  means  of  the  above  process;  how¬ 
ever,  the  sample  space  now  spans  the  discrete  time  interval  from  time  k-n  to 
time  k+1.  The  sample  set  X  of  equation  (53)  is  now  defined  as 


A  new  weighted-least-squares  solution  for  the  straight-line  regression  fit 
through  the  n  samples  is  found  and  the  process  repeats  as  outlined  above • 

Note  that  the  solution  of  the  nonlinear  relationships  of  (54),  (57),  and 
(67)  are  obtained  from  an  iterative  procedure.  Equations  (54),  (57),  and  (67) 
are  nonlinear  as  a  result  of  the  bisquare  weight  function,  W ^ ,  given  by  the 
relationship  of  equation  (60). 

SIMULATION  RESULTS 

The  system  of  figure  4  (with  a  -  2.0  and  b  -  3.0)  was  simulated  on  a 
digital  computer;  a  simulated  sensor  monitored  the  output,  x.,(t),  where  the 
output  was  measured  discretely  in  time  and  corrupted  by  sensor  noise,  v(k). 

The  measurement  noise,  v(k),  was  zero  mean  with  a  variance  of  25.  The  system 
was  driven  by  a  deterministic  forcing  function  u(t).  No  process  noise,  w(t), 
was  added  to  the  system  forcing  function.  The  deterministic  forcing  function 
was  a  pulse  with  a  duration  of  22  seconds  and  a  magnitude  of  500  units,  as 
illustrated  in  figure  5.  Also  shown  in  figure  5  are  records  of  the  values,  as 
functions  of  time,  of  the  state  variables  x2(t)  and  x^t).  The  output  x^t) 
was  sampled  at  a  rate  of  five  times  per  second. 


TIME.  Mcortch 


Figure  5.  System  forcing  function  and  values  of  state  variables. 


A  conventional  Kalman  filter,  without  any  a  priori  knowledge  of  the 
forcing  function,  u(t),  or  the  time  at  which  the  forcing  function  was  initi¬ 
ated,  was  used  to  process  the  measurement  data  z( k) •  The  conventional  Kalman 
filter  did  not  detect  the  influence  of  the  deterministic  forcing  function  on 
the  state  variables,  as  illustrated  in  figure  6.  Since  a  priori  data  dictated 
that  there  was  no  process  noise,  the  elements  of  the  Kalman  filter  gain  matrix 
associated  with  the  observed  variables  approach  zero;  thus  the  estimation 
process  has  severed  itself  from  the  measurement  process  and  ignores  new  data 
brought  forth  by  additional  measurements. 

ESTIMATES  OF  VALUES 
OF  STATE  VARIABLE 


10  12  14  16  18  20  22  24  26  |  28  30  32 

TIME,  seconds 


Figure  6.  Estimation  using  adaptive  Kalman  filter  without  robust 
smoothing  compared  to  estimation  using  nonadaptive  conventional 
Kalman  filter. 

When  the  elements  of  the  Kalman  filter  gain  and  covariance  matrices  are 
weighted  by  the  adaptive  procedure  outlined  above  (sample  statistics  of  the 
innovations  sequence  are  used  to  adapt  the  respective  weights),  the  filter  no 
longer  divorces  itself  from  the  measurement  process.  Additional  data  brought 
forth  by  the  measurement  process  are  used  to  update  the  estimates  of  the  state 
variables.  This  is  illustrated  in  figure  6.  However,  since  the  adaptive 
procedure  uses  sample  statistics,  the  estimates  contain  periodic  outliers. 


The  filter  will  run  for  a  period  of  time,  then  monitor  the  innovations  se¬ 
quence  to  update  the  adaptive  weights.  It  is  this  monitoring  of  the  innova¬ 
tions  sequence  to  obtain  new  information  which  causes  the  periodic  outlier  to 
appear  in  the  estimates.  The  subsequent  processing  of  the  adaptive  estimates 
by  a  robust  smoother  reduces  the  level  of  mean  square  error  and  the  periodic 
outliers . 

The  smoothed  estimates  of  the  measured  output  state  variable,  x^(k),  are 
shown  in  figure  7.  Figure  8  presents  an  overlay  of  the  records  of  figures  6 
and  7  over  a  5-second  expanded  time  interval. 

CONCLUSIONS 

State  variable  estimation  in  the  presence  of  unknown  a  priori  system 
information  (noise  statistics,  forcing  functions,  and  system  dynamics)  is  not 
an  easy  problem.  There  are  no  clear-cut  solutions.  This  paper  addresses  only 
one  of  the  above  problems  (no  information  about  the  system  deterministic 
forcing  functions).  The  concepts  presented  relative  to  this  particular  prob¬ 
lem  address  the  limited  class  of  linear  system  dynamics  with  associated  linear 
measurements.  Nonlinear  system  dynamics  with  associated  linear  measurements, 
however,  are  not  precluded. 

Intimates  of  the  state  variables  using  the  adaptive  process  for  the 
system  during  the  periods  when  the  system  is  not  being  forced  are  relatively 
close  to  those  of  the  conventional  Kalman  filter  for  congruent  periods,  but 

there  is  some  degeneration  because  the  estimator  is  no  longer  optimal.  During 

the  periods  when  the  system  is  being  forced,  a  vast  improvement,  as  compared 
with  those  estimates  of  the  conventional  Kalman  filter,  is  realized  with  the 
adaptive  gain,  covariance  weight,  and  associated  robust  smoothing  procedure. 
The  estimates  derived  with  the  adaptive  procedure  during  the  periods  of  system 
forcing  do,  however,  contain  a  considerable  level  of  mean-square  error.  This 
seems  to  be  a  prevailing  shortfall  of  adaptive  estimation  procedures.  The 

tradeoff  is  knowing  more  about  the  values  of  the  state  variables  (less  mean 

error)  against  more  mean  square  error  in  the  respective  estimates. 


21 


6  8  10  12  14  16  18  20  22  24  26  28  '  30  TIME,  seconds"" 

Figure  7.  Estimation  using  adaptive  Kalman  filter  with  robust  smoothing. 
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Figure  8.  Comparison  of  filtering  techniques  with  and  without  robust  smoothing. 
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APPENDIX  A 


DERIVATION  OF  RECURSIVE  ESTIMATORS  FOR  SAMPLE  STATISTICS 


Presented  in  this  appendix,  the  derivations  of  the  recursive  estimators 
for  the  sample  mean  and  sample  variance  based  on  N^  observations. 


RECURSIVE  SAMPLE  MEAN 


The  expression  for  the  recursive  estimator  for  the  sample  mean  at  time 

t  is 
k-1 


-  i\  i 

n(k  "  ^  2  n(j) 

l  j=k-N 


(  A .  1  ) 


The  expression  for  the  recursive  estimator  for  the  sample  mean  at  time  t 


~~  2  n(j] 

1  j*k-N  >1 


(A. 2) 


Equation  (A.1)  is  subtracted  from  (A. 2)  to  give 


n(  k)  -  n( k  - 


r  k  k-1  -1 

D  *  7T-  2  n<  -  2  n(j)  ' 

t  L  j=»k-Nt*1  j*k-Nt  J 


(A. 3) 


When  the  terms  under  the  sumnation  of  (A. 3)  are  expanded  out  and  appropriate 
cancellations  take  place,  the  recursive  estimator  for  the  sample  mean  at  time 
t  is  given  as 


n(k)  «  n(k  -  1 )  ♦  “  [n(k)  -  n(k  -  N.  )  ] 

Nl  1 


(A. 4) 
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RECURSIVE  SAMPLE  VARIANCE 


V*‘»  S  v" 


Expanding  the  terms  under  the  summation  gives 


(A.  11) 


